Fluctuations in fluid invasion into disordered media 
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Interfaces moving in a disordered medium exhibit stochastic velocity fluctuations obeying universal 
scaling relations related to the presence or absence of conservation laws. For fluid invasion of porous 
media, we show that the fluctuations of the velocity are governed by a geometry-dependent length 
scale arising from fluid conservation. This result is compared to the statistics resulting from a 
non-equilibrium (depinning) transition between a moving interface and a stationary, pinned one. 
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The dynamics of interfaces, domain-walls and fronts 
in disordered systems has attracted tremendous attention 
not only because of practical applications but also for the 
general interest in the emerging equations of motion and 
scaling laws. The presence of randomness is felt through 
both frozen (or "quenched" ) noise and "thermal" fluctua- 
tions. Experimentally tested scenarios include fire fronts 
in combustible media [l[ , cracks propagating in solids 
@- 0, 0] , domain walls in magnets [S IS 0> IS La] > an d in- 
terfaces in multi-phase flows ( "imbibition" or "capillary 
rise" phenomena) [IS EL IS IS 14 1, to name a few. 

In these examples the fluctuations of the fronts 
can generally be described via self-affine fractals, with 
dynamics exhibiting critical correlations. For time- 
independent randomness, such as the pore structure in 
fluid invasion, it is possible to concentrate on the zero- 
temperature phase-transition which separates moving in- 
terfaces from pinned ones. The situation is then repre- 
sented in terms of an order parameter, the velocity of 
the interface, and a control parameter, the driving force. 
The interface is described by a height h(x.,t), at posi- 
tion x € R d and time t, with its dynamics given by a 
Langevin equation dth = £({h},a({h},x,t)). The kernel 
C contains the deterministic contribution of the interface 
and the random contribution of the noise configuration 
a({h}, x, t). Both may depend non-locally on the entire 
interface configuration, denoted by {h}. 

The scaling exponents of /i(x, t) are often studied but 
are by no means the only interesting feature of this prob- 
lem. It is also possible to study the amplitude of the 
various correlation functions, together with the persis- 
tence properties of the interface [lq . or the scaling of 
probability distributions of quantities such as deviations 
of the interface from its average position or velocity. 

Here, we study the velocity fluctuations of driven fluid 
interfaces in disordered porous media in an imbibition 
situation, where a viscous fluid displaces air or a less vis- 
cous fluid. It is an important problem in the general 
field of fluid flow in porous media, with many engineer- 
ing and technological applications. Fluid invasion and 
multi-phase flows are often analyzed via lattice models 
of pore networks (e.g. for oil industry applications), or 



via "upscaling techniques" , which summarize the physics 
on a coarse-grained (volume element) scale. One often 
studied aspect is the averaging of the scale-dependent 
permeability k(t), from the smallest pore scale up to ge- 
ological scales IS, II , IS , IS] . It is not trivial to arrive to 
a coarse-grained description of these systems, since the 
physics is governed by the (conserved) fluid flow in the 
porous bulk coupled to microscopic interfaces, menisci, 
between the two fluids. 

We report on the behavior of the spatially averaged in- 
terface velocity v(t), directly related to the liquid intake 
under an externally controlled pressure. We first show 
that there exists a scaling relation between the veloc- 
ity fluctuations Av and the average velocity v — (v(t)). 
This scaling is also reflected in the power spectrum of 
v(t) as we discuss in detail. The scaling arguments are 
backed with numerical simulations of a phase-field model 
for imbibition 2fj. We also study the interface velocity 
fluctuations for non-conserved, "local" dynamics, where 
we find a different scaling relation , ag ain backed by sim- 
ulations of the appropriate model [21| . Both the average 
flux as well as its fluctuations are experimentally well ac- 
cessible so our theoretical approach allows for an easy 
assessment of the type of dynamics. The detailed exper- 
iments of the Barcelona group and others point into that 
direction [liillll EE! fill ] , and the main result extends 
to other systems where the front invasion depends on a 
conserved quantity. 

The interface advances because fluid is transported 
from the reservoir through the medium. By Darcy's law, 
fluid flow is proportional to a pressure gradient, 



j = - VP. 

V 



(1) 



Here r\ is the viscosity of the liquid and k the permeabil- 
ity of the medium. Deviations from Eq. ([1]) are known to 
arise for microscopic reasons, such as, e.g., inertial effects 
when a meniscus enters a pore, clogging by advected par- 
ticles as in the case of a dye, and the partial wetting of 
capillaries. Nevertheless, the linear form already leads to 
rich phenomena [22j |. 

For an incompressible fluid, a Laplace-equation 
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V 2 P = governs the pressure field in the bulk. The 
appropriate boundary conditions are (i) P = Pr at the 
contact to a reservoir of liquid with constant pressure, 
and (ii) Pi = j*K,+P c +Pq at the interface. Condition (ii) 
is a superposition of atmospheric (pj) and capillary (P c ) 
pressure on the pore scale and the effect of curvature K, 
on a coarse-grained scale by an effective interface tension 
7*. For weak disorder, as arise in e.g. paper, it is close 
to the bare value 7, but may be much smaller in porous 
materials with weakly connected channels such as Vycor 
glass [231 ] . 

For Pr = Po imbibition is spontaneous with an average 
velocity v = d t h— (n/ij) Pc/h, where h(t) is the average 
height. This leads to Washburn's law of capillary rise, 
h^t 1 / 2 24 1 . Forced imbibition arises by increasing the 
pressure at the reservoir as the interface advances. A 
constant mean pressure gradient, |VP| = (Pr — P\)/h = 
tjv/k, keeps the average velocity v constant. 

Quenched disorder in the porous structure roughens 
the interface around its average position. Capillary dis- 
order p c (r) = P c + Sp c (r) acts only at the interface, and 
permeability disorder, n(r) = k + 5k(t) controls the flux 
of liquid from the reservoir to the interface. If ro is the 
typical pore size and Sr^ its deviation, then p c ~ 7* /ro 
and K^r 2 ., with 5p c / 5n^p c / k [22I [25|. 

Effective (capillary) interface tension and the average 
pressure gradient smoothen the interface and allow for 
correlated roughness only up to a lateral length scale 



qv 



(2) 



related to the capillary number Ca = rjv/j* 2^, 22|. Its 
presence and dependence on v has also been confirmed 
experimentally [lfj, |l3j . The ratio of disorder strengths 
imposes another length scale 



vrj Sk Ca 



(3) 



separating regimes of different roughness. On scales 
shorter than £ K interface roughness is caused by capil- 
lary disorder, on larger scales by permeability disorder 
[22L [25 1 . We consider the interesting case of a slowly 
advancing front, where £ K > £ c , and capillarity induced 
fluctuations prevail, and one has a well defined interface 
(for stability a Ca below ~ 0(1) is needed (26[). Then, 
pcrmeablity noise and its possible correlations with cap- 
illary disorder can be ignored. 

The scaling behavior in imbibition can be demon- 
strated by model simulations [2fJ. A scalar phase field 
cj)(r, t) denotes by <fi = 1 the invaded (wet) and —I the 
non-invaded (dry) regions respectively. An energy func- 
tional 



m = J < 
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couples the phase field to a space dependent quenched 
randomness a(r) with average a and standard deviation 
Aa. This term models capillarity disorder, with sign 
(a > 0) favoring the wet phase 0=1. The model is defined 
on the half-space, r = (x, y) with y>0, with boundary 
conditions (/)= 1 at y — (coupled to a wet reservoir at 
the bottom) and initial condition (f>(r, t — 0) = — 1 (dry). 
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FIG. 1: Left: Interfaces h(x,t) in a system of width L = 
500 and average separation h = 320 from the reservoir. The 
configurations are separated by time intervals At — 500 at 
240 different times. Interface propagation is made visible by 
adding the constant shift velocity v. The sequence of 20 black 
interfaces shows pinned and moving regions. Right: time 
series of v(t), the velocity for L = h = 40, 80, 160 .. . 1280 
from top to bottom, logarithmic scale in v(t). 
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2: Activity pattern in space (vertical) and time (horizon- 

L = 256, h — 160, lateral periodic boundary conditions. 
High velocities dark, low velocities light, nonlinear grey-scale 
for better visibility. At any given time activity is restricted 
to n = L/h ~ 2 narrow avalanches of width £ c . 

The corresponding chemical potential \i = b~T /5(j) (with 
the role of pressure in the context of imbibition) drives a 
current j = — kV/i with a rescaled and possibly spatially 
varying mobility k(r). By local conservation of this 
yields the dynamical equation 

d t 4> = -V • k(r) V [V 2 </> + <j) - (j) 3 + a(r)] . (5) 

With this model, a Young-Laplace relation between the 
chemical potential at the interface and the curvature, 
Mint —IC — a arises naturally. 

The numerical simulations are performed, in a way 
that mimics forced-flow imbibition, by continuously 
shifting the fields (j>(r,t), R(r) and a(r) downward at ve- 
locity v, as in the experiments of Ref. [27|. The station- 
ary interface then has average height h — a/(2v). For 
the Ca at hand the permeability disorder has no influ- 
ence (see [22}), nor have correlations between k(r) and 
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a(r). Only simulations with capillary disorder are thus 
discussed here. Fig.[T]shows examples of model interfaces 
in d+1 = 2 dimensions, together with timeseries of average 
velocity or liquid intake, v(t) exhibiting typical fluctua- 
tions. Fig. [3 shows a spatiotemporal activity pattern of 
the local interface velocities d t h(x,t). 

The dynamics of the roughening interface can be un- 
derstood as a sequence of avalanches. At any given time 
only a small localized portion of the interface is mov- 
ing fast. There the Laplacian pressure field drags liquid 
from a surrounding region of lateral size h d , region which 
therefore remains pinned. Further away the interface is 
not affected and moves independently. Within this re- 
gion, the moving part has maximal capillary forces p c . 
However, it eventually encounters disorder with lower p c 
and gets pinned. At this point, another region elsewhere 
starts to move. Propagation by avalanches in imbibition 
and their size distribution with one fixed scale has been 
found and discussed in experiments 14j. It is visible in 



Figs. Q] and [H The size of these moving avalanches is 
given by the cutoff lengthscale £ = £ c or £ K , The regime 
governed by capillary disorder has global roughness ex- 
ponents x~ 1-25 and a different local exponent xioc = 1 in 
d+1 = 2 dimensions [io|, [H, [Hj] . For higher Ca than those 
considered here, permeability disorder starts to play a 
role and the exponents are x = Xioc = 1 in d+1 = 2 dimen- 
sions UlSil. 

A key feature is the relation between an avalanche 
duration r and its volume s(r). In general, a scaling 
s(t) ~ r 7 , may be expected, with exponent 7 determined 
as follows. Consider a region of lateral size I swept over 
by an avalanche. Its vertical extent is related to I by the 
local roughness exponent, w~£ Xlaa [13, [Hj], with volume 
s ~ £ d +xioa _ ^he avalanche is driven by a higher capil- 
lary pressure in the moving region as compared to other 
parts of the front. Since locally the values of p c are inde- 
pendent random quantities, the excess driving force (the 
difference of typical p c and the large fluctuation) and ve- 
locity decrease with size as v ~ £~ d / 2 . Therefore, the 
avalanche sweeping time scales as t = w/v <~ pt^+ d / 2 
and leads to the relation s ~ r 7 with 



Xloc 



Xioc + d/2 



(G) 



For d — 1 we have xi oc = 1 [20( and expect 7 = 4/3. In 
rf = 2 wc would expect using xi oc = X =3/4 [29| 7 = 11/7. 

To see the scaling of velocity fluctuations Av we first 
relate the average avalanche velocity w ava to the overall 
velocity v by the relative size of moving parts, w ava ~ 
(h/^) d v. Next, we assume a simple relation Au ava ~u ava 
between the average avalanche velocity and its fluctua- 
tions, justified by the slow, intermittent motion. Thus 
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by factors resulting from independence on scales larger 
than h and the relative fraction of moving parts on 
smaller scales. The relation between average avalanche 
velocity and fluctuations yields in general, both for spon- 
taneous and forced imbibition 



Av 



(8) 



In spontaneous imbibition with rising height v~l/h and 
Av^v 1 ~ d / 2 L- d / 2 . Eq. © thus predicts that the fluctu- 
ations depend on the geometry (via h and L) and on the 
average velocity. 

In the left panel of Fig. [3] we show the distributions 
of P(v(t)) obtained in phase field simulations for several 
values of v, and constant square aspect ratio L = h: They 
can be scaled by an ansatz P(v) = \P{^-^)- The shape 
of P(v(t)) resembles a Gumbel distribution of extreme 
value statistics reflecting avalanches in regions of maxi- 
mal capillary pressure, but our data are too scarce for a 
careful analysis. 
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FIG. 3: Left: Rescaled distributions -p(<£±) f or 

various 

values of l/voch = L = AQ, . . . , 1280. Right: Power spectrum 
S(w) = (\v{u)\ 2 ) oc w" 4/3 for h = L= 1280. Same exponent 
relates average avalanche size s and duration r. 

The avalanche motion is reflected in the power spec- 
trum of the interface velocity, S(v) — (\v(uj)\ 2 ) ~ oj~ a . 
The exponent a may equal the exponent 7 discussed pre- 
viously, if the avalanches have a self-affine fractal spa- 
tiotemporal structure [131 . The right panel of Fig. [3] com- 
pares the power spectrum S(u) = (\v(u))\ 2 ) of the mean 
interface velocity (solid line) to the average size s(r) of 
avalanches with duration r (dashed line). Both follow 
power laws, uj~ a and r 7 with a = 7 « 4/3, as predicted 
by Eq. ©. 

The left panel of Figure 2] shows Aw vs. v for two dif- 
ferent choices of the system geometry: h varying, with L 
fixed, and with L~h. We observe Aw oc v a with a = 1 for 
a square shape and a= 1/2 for fixed width in agreement 
with our scaling picture. We checked also that for L^>h 
the distribution P(v(t)) gets sharper with increasing L, 
in agreement with the implication of the central limit 
theorem for independent random variables Av ~ L^ 1 ! 2 
(Eq. ©). 

An interface driven by a force F in a random medium 
without any conservation law obeys different scaling. In 



4 



the presence of a critical (zero-temperature, neglecting 
thermal creep) value F c separating pinned and unpinned 
propagation a correlation length £~AF _l/ ensues where 
AF = F — F c and v is the correlation length exponent. 
Meanwhile, the order parameter follows the generic scal- 
ing v~ AF . In the critical region, when AF is so small 
that £ « L, the order parameter exponent 9 takes a non- 
trivial value smaller than unity. For slightly larger driv- 
ing forces, we have L > £. In the simulations here we 
considered L/£ = O(10), so we have v ~ AF or a trivial 
effective 9 e g = 1 . 

In a system of lateral size L there are N = (L/^) d 
independent sub-volumes. Inside these, the fluctuating 
locally averaged velocities v\ oc are independent random 
variables with mean v and variance Sv 2 . Their average, 
the instantaneous velocity v(t) = 'Yl li v l ° c /N then has 
fluctuations measured by 



Av ~ Sv/Vn ~ v 1 -^/^ /Via. 



(9) 



Numerical simulations are performed with a cellular 
automaton for the QEW equation [21|, dh/dt — W 2 h+ 
F + T)(x,h(x,t)), where r] is a quenched noise term with 
a bare correlator (-q(r)ri(r')) =D6(r — r') with strength 
D. r is a surface tension similar to 7*. The data are 
presented in the right panel of Fig. |4j They fit reasonably 
to the scaling ansatz of Eq. with 1— dv/(26 e g) = l/3, 
derived from the known critical value z/«4/3 of the one- 



dimensional depinning transition 21| and 9 c g = 1 
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FIG. 4: Left: Fluctuations Av vs. mean v for constant width 
£ = 500 (open circles) and L = h (full circles), follows Av^v a 
with a = l/2 and 1 resp., cf. Eq. (JSJ. Lines with slope 1/2 and 
1 resp. guide the eye. Right: Same for local model of [2lj . 
with system size L — 10 4 . Line has slope 1/3. 

To summarize, we have presented a scaling theory for 
interface velocity fluctuations. Without detailed mea- 
surements of interface configurations but by the easily 
accessible fluctuations of the average interface propaga- 
tion one obtains information about the universality class 
of the process. Our arguments show the difference be- 
tween dynamics with and without a bulk conservation 
law behind the interface. This is also manifest in the 
power spectrum of v(t), where the avalanche-like dynam- 
ics allows to understand the ensuing 1 //-noise. The main 
argument applies to any dimensionality and is not depen- 
dent on the detailed model. A similar reasoning should 
extend to other systems where the global conservation of 



a quantity is important. One experimental possibility is 
vortex penetration dynamics in (high-T c ) superconduc- 
tors, — for which there is no known coarse-grained de- 



scription 



31, 3 



since the vortex density is conserved. 
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